Localization and Kosterlitz-Thouless Transition in Disordered Graphene 
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We investigate disordered graphene with strong long-range impurities. Contrary to the common 
belief that delocalization should persist in such a system against any disorder, as the system is ex- 
pected to be equivalent to a disordered two-dimensional Dirac Fermionic system, we find that states 
near the Dirac points are localized for sufficiently strong disorder and the transition between the 
localized and delocalized states is of Kosterlitz-Thouless type. Our results show that the transition 
originates from bounding and unbounding of local current vortices. 
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It is well-known that the electronic spectrum of 
graphene can be approximately described by relativis- 
tic Dirac Femions[ll, This is due to the linear dis- 
persion relation at low energies near two valleys associ- 
ated to two inequivalent points K and K' at the cor- 
ner of the Brillouin zone(3|- The relativistic dispersion 
gives rise to several remarkable phenomena. Unlike non- 
relativistic Schrodinger fermions in two dimensions 
Dirac fermions cannot be trapped by a barrier due to 
the Klein paradox, a property of relativistic quantum 
mechanics 4]. Theories based on the two-dimensional 
(2D) single flavor Dirac Hamiltonian also predict that 
Dirac fermions cannot be localized by disorder 

The great majority of experimental and theoretical 
studies of graphene 0, has focused on the effect of the 
relativistic electronic dispersion on different phenomena 
such as Landau level structure or quantum Hall ferro- 
magnetism. However, the validity of single flavor Dirac 
fermion picture for disordered graphene is only approxi- 
mate and relies on two premisses: (1) The spatial range 
of the impurities is long enough to avoid inter-valley 
scattering[5| (for short-range impurities, strong inter- 
valley scattering can lead to localization^, [ll|, [IjI); (2) 
Even in a disordered graphene with long-range impurities 
that completely suppress inter-valley scattering, a single 
valley Dirac Hamiltonian is only valid when weak impu- 
rities are considered. Since the approximated linear rela- 
tivistic dispersion is valid near the Dirac valley points, K 
and K', the approximation cannot be carried out in a re- 
gion with a strong impurity where the deviation from the 
Dirac point is large enough so that higher order correc- 
tions to the energy spectrum become relevant 'B'l. There- 
fore, the application of single valley Dirac Hamiltonian 
to disordered graphene is limited to weak long-range im- 
purities. Indeed, localized states in disordered graphene 
near Dirac points have been observed experimentally [13] 
and numerically [3, All the above beg the physical 
question: how does graphene behave in the presence of 



strong long-range impurities? 

In this letter we investigate several novel phenomena 
induced by disorder with strong long-range impurities in 
graphene. We calculate the scaling properties of dis- 
ordered graphene in the framework of a tight-binding 
model and finite-size scaling. Instead of delocalization 
we find that, in the presence of strong long-range impu- 
rities, states near the Dirac points are localized. Local- 
ization arises from enhanced backscattering due to the 
deviation from linear dispersion in the strong impurity 
regime. We show that there is a metal-insulator tran- 
sition (MIT) as a function of the disorder strength and 
chemical potential. On the delocalized (metallic) side, 
the conductance is independent of the system size, which 
is a characteristic of the Kosterlitz-Thouless (K-T) [3] 
type transition in conventional 2D systems with random 
magnetic field 18, l^l or correlated disorder [^^j- We ver- 
ify the Kosterlitz-Thouless transition nature of the MIT 
by explicitly identifying the bounding and unbounding 
vortex-anti-vortex local currents in the system. 

The TT electrons in graphene are described by the tight 
binding Hamiltonian (TBH) 



i7 = ^y.4c, + i^(4c,-fH.c) 



(1) 



where c| (ci) creates (annihilates) an electron on site i 
with coordinate r^, t (~2.7eV) is the hopping integral be- 
tween the nearest neighbor carbon atoms with distance 
a/^/3 {a ~ 2.46A is the lattice constant), and Vi is the 
potential energy. In the presence of disorder, Vi is the 
sum of contributions from Nj impurities randomly cen- 
tered at {r^} among N sites Vi = Ylm=i exp(-|ri - 
rmp/(2^)), where Um is randomly distributed within 
(—14^/2,14^/2) in units of t. Different random configu- 
rations of graphene samples with same size, ^, W and 
n.i = Nj /N constitute an ensemble with definite disorder 
strength. This model has been widely used in investigat- 
ing the transport properties in graphene 21 , 2^ 3] ■ 
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At zero temperature, the two terrainal dimensionless 
conductance gL of the sample between perfect leads at 
Fermi energy Ep can be written in terms of Landauer- 
Biittikcr formula [Hi: 



(2) 



where t is the transmission matrix and the factor 2 ac- 
counts for spin degeneracy. Equation ([2]) can be numer- 
ically evaluated by recursive Green's function method 
[ist for systems with rather large size. For the pur- 
pose of scaling, the contact effect should be subtracted 
from gL to yield the "intrinsic conductance", g, defined 
as l/g — l/gL — l/(2iV(7), where Nc is the number of 
propagating channels at Fermi energy Ep and l/(2A^c') 
is the contact resistance [i^. The conductance g then 
receives contributions solely from the bulk and thus has 
the same scaling property as if it were obtained by the 
transfer matrix method [l^]. The scaling function (ol. 



/3 = 



d(ln.g) 
dlnL 



(3) 



(. . .) being the average over random ensemble, is used to 
determine the localization properties; /? < and /3 > 
correspond to the insulator and the metal, respectively. 

We plot the size dependence of (lug) with ^ — 1.73a, 
m = 1%, for different Ep and W in Fig. [1] The sam- 
ples are set to be square shaped with length L. Peri- 
odic boundary conditions in the transverse direction are 
adopted to exclude the edge states of the zigzag edges p2|. 
The potential range ^ here is chosen long enough to avoid 
obvious inter- valley scattering 0, [2§|, and the scaling 
S,/L^ - is irrelevant. When Ep < E^ = O.lt (Fig. 
□ (a)) or > Wc = 2t (Fig. [T] (b)), (g) is monotonically 
decreasing with increasing L, which means the wavefunc- 
tions are localized. Otherwise, when Ep > Ec = O.lt 
(Fig. [l](a)) orW <Wc = 2t (Fig. [T](b)), (Ing) curves for 
different sizes merge, suggesting a delocalized state with 
finite conductance in the thermodynamic limit. However, 
they are not real metals with (3 > 0. All the states with 
W G (0, Wc) are within the metal-insulator transition 
(MIT) region with (3 = 0. Even in the cases of extremely 
weak disorder with W = 0.25t and W = O.lt (see the 
inset of Fig. [D (b)), except for a vanishing even-odd like 
fluctuation, {\ng){L) doesn't seem to show a tendency 
to be increasing nor decreasing. In Fig. [21 the universal 
/3(ln g) is plotted from the same data in Fig. [TJ showing a 
critical conductance In ~ 1 separating the delocalized 
states with /? = and localized states with f3 < 0. This 
phenomenon corresponds to a disorder-driven Kosterlitz- 
Thouless (K-T) type transition that has been observed 
in many disordered 2D systems 17, IH, [H, As can 
be seen from Fig. [1] (a) and the phase diagram (inset of 
Fig. [2]), states in the low energy region are more easily 
localized. 

The existence of localized states near the Dirac point 
is in contrast to the belief that Dirac fermions are robust 




FIG. 1: (Color online) The scaling of conductance for long 
range disorder (^ = 1.73a, n/ — 1%): {a.){\ng) as func- 
tions of the Fermi energy Ep with fixed disorder strength 
W = 2t(note: the bandwidth is 6t)); (b) (Ing) as functions of 
disorder strength W with fixed Fermi energy Ep = O.lt. The 
insets are the same data plotted as functions of size L. Each 
(Ing) is an average over 100 ~ 400 random realizations. 



against localization, especially in the presence of long 
range impurities that can effectively prohibit the inter- 
valley scattering. In order to gain insight in the nature 
of the localization transition, we now turn back to the 
dispersion structure of realistic graphene. In the absence 
of disorder {Vi = 0), the upper (+) and the lower (— ) 
bands touch at two Dirac points K = (|^, ) and 

K' = (if , --j^)- When \E\<t, the dispersion consists 
of two valleys centered at Dirac points. Near each Dirac 
point, e.g. K, the energy bands can be expanded as[^ 

i?±(q) = ±^|q| ± ^sin(3a(q))|qp + 0{q% (4) 

where q = k — K is the momentum measured from K 
and a(q) € [0, 27r) is the angle of vector q. The first 
term in the r.h.s of ^ corresponds to the Dirac Hamil- 
tonian, but non-linear terms will be prominent when q 
(or E) is increased. Even when \E\ < t, the second term 
in (|4]) can lead to non-trivial consequences. First, the 
quadratic dependence on momentum |qp gives rise to 
non- vanishing backscattering probability and thus a ten- 
dency to localization. Second, the angular dependence 
factor sin (3Q;(q)) ("trigonal warping") breaks the perfect 
symmetry of the cone-like valley. The pseudo-time rever- 
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FIG. 2: (Color online) The scaling function 13 — ^i^j^ 
tained from the data in Fig. [T] The inset is the schematic 
phase diagram for ^ = 1.73a, ni = 1%. 



sal symmetry [s, 22| restricted to each valley is destroyed. 
When \E\ > i, the linear approximation and double- 
valley structure collapses completely. Although graphene 
cannot be experimentally doped to a bulk Fermi en- 
ergy far away from the neutral point (Dirac points) (e.g., 
Ep ^ t), the local potential of impurities might still be 
high enough to create non-Dirac scatterings. The ob- 
served localization originates from the non-Dirac behav- 
ior due to higher order corrections to the dispersion. 

To confirm this, let us consider the simplest case of a 
single long range impurity in the center of a graphene 
sheet. If inter- valley interaction is effectively prohib- 
ited and the regime of Klein tunneling[4] holds, there 
should be no bound states, no matter how high the po- 
tential barrier is. After diagonalizing the Hamiltonian 
for a graphene sheet with TV sites, the spatial extension 
of eigenstate \tpn) = Y^iLi o-nicl\0) with eigen-energy En 
can be characterized by the participation ratio 



N 



N 



^n = (E«-)V(^E«™)' 



(5) 



which is a measure of the portion of the space where the 
amplitude of wavefunction differs markedly from zero. 
For an extended state, R has a finite value (typically 
close to 1/3 in the presence of disorder), whereas for a 
localized state R approaches zero proportional to {l/N) 
[2^ . The results for ^ = 1.73a with different potential 
height V >Ois plotted in Fig. O For small V (Fig. [3] 
(a) and (b)), where the electronic behaviors inside and 
outside the barrier are Dirac-like (Fig. [3] (e) and (f)), 
there are no bound states. When V is increased, bound 
states with small R begin to appear in the negative en- 
ergy region near the Dirac point, as seen in Fig. [3] (c). 
For positive injected energy (orange arrow with solid line 
in Fig. [3] (g)), the electron is not far from K both in- 
side and outside the barrier, so the regime of Klein tun- 
neling is still valid and the electron cannot be trapped. 




FIG. 3: (Color online) Left column ((a)^(d)): The partici- 
pation ratio R as functions of energy E for a graphene with 
= 70 X 40, in the presence of a single impurity at the 
center with ^ — 1.73a (long range) with different potential 
height V > 0. Right column ((e)— >(h)): Schematic diagrams 
of scattering process corresponding to their left counterparts. 
The electron is injected from the left, scattered by the barrier 
at the center and eventually transmitted to the right (thick 
orange arrows). The dispersion configurations in these three 
regions around K are plotted, where the red lines mark the 
ideal Dirac dispersion E±(k) = ±^|k-K| and the olive part 
is that for graphene calculated from TBH. Discrepancies be- 
tween them at high energy can be clearly seen. 



On the other hand, for negative energy (orange arrow 
with dashed fine in Fig. [3] (g)), the electron sees a non- 
Dirac barrier (pointed by the yellow arrow) . This causes 
strong back-scattering and localization around the impu- 
rity. When V is increased further (Fig. [3](h)), even elec- 
trons in the positive Dirac region will encounter strong 
back-scattering in the barrier and will be localized (Fig. [3] 
(d)). For negative V (not shown here), all the results are 
similar, except that the localized states now first appear 
in the positive region near Dirac point. In conclusion, 
the localized states originate from back-scattering at the 
barrier due to its deviation from Dirac behavior at the 
Fermi level. The states near the Dirac point with low 
density of states will be more sensitive to backscattering 
and will be localized fi rst, as in the case of conventional 
disordered systems 25l |26|, [27 1 . 

Why is the MIT in disordered graphene of K-T type? 
The K-T transition is a typical topological transition 



which has been understood as unbounding of vortex-anti- 
vortex pairs [l6l |. For instance, in the high temperature 
phase 2D XY model, a plasma of unbounded vortices and 
anti- vortices of local spins gives rise to an exponential de- 
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anti-vortices (n < 0) are closely bounded, corresponding 
to the "low temperature" phase of 2D XY model with 
quasi-long range correlations. In the localized phase (Fig. 
|4] (b)), there are a large number of current vortices and 
anti- vortices. Many of them are unbounded, correspond- 
ing to the "high temperature" phase of 2D XY model 
without long range correlations. This offers an explicit 
picture of the microscopic origin of the disorder driven 
K-T transition in graphene. 

In conclusion, we find a Kosterlitz-Thouless type 
metal-to-insulator transition as a function of disorder 
strength or Fermi energy in disordered graphene with 
strong long-range impurities. We explicitly demonstrate 
the KT nature of transition by showing the bounding 
and unbounding of local current vortexes. One unique 
feature about the K-T transition is a scaling of expo- 
nential form near the transition point. The conductance 
g oc exp(— ■^^^==) where a is a constant. Recently, 
the MIT near the neutral point of graphene has been 
observed in graphene nanoribbons[l3j which are quasi- 
one-dimensional systems. Our results can be tested in 
experiments with large nanoribbon radius. 

We thank D.-X. Yao, W.-F. Tsai and C. Fang for useful 
discussions. YYZ and JPH were supported by the NSF 
under grant No. PHY-0603759. 



FIG. 4: (Color online) Typical configurations of local currents 
in (red arrows) and potential Vn (color contour) on two sides 
of K-T type MIT with iV = 56 x 32 sites, = 1.73a, m = 1% 
and Ef = O.lt. (a): W = l.lt (delocalized); (b): W = 2M 
(localized). The size of arrows is proportional to the logarithm 
of current value. Carbon hexagons with topological charge 
n ^ are marked explicitly with blue numbers. Both plots are 
in the same random realization of impurities, with diflterent 
potential height therefore eflPectively different W. 



cay of spin correlation function; in the low temperature 
phase, vortices and anti- vortices arc bound to each other, 
leading to a power law correlation function. This can be 
clearly seen in the present problem if the local currents 
are identified with local spins in XY model. 

The bond current vector ii^rniEp) per unit energy 
pointing along the bond between sites I and m can be 
calculated using Green's functions 0, [1^ H^]- It is 
more convenient to investigate the "current flow vector" 
i; = h^m defined on site I, where the vectorial sum- 
mation is taken over the nearest neighbors of site I [2^ . 
The current flow i; is a vector with angle di G [0, 27r). The 
topological charge n of local currents on a closed path 
can now be defined as usual: n — § \79 ■ d£. In Fig. 
m typical distributions of local currents on both sides of 
MIT are plotted. As expected from the K-T picture, in 
the delocalized phase (Fig. H] (a)), vortices (n > 0) and 
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